Vis Comput 
DOI 10.1007/s0037 1-01 1-0610-y 


ORIGINAL ARTICLE 


Harris 3D: a robust extension of the Harris operator for interest 


point detection on 3D meshes 


Ivan Sipiran - Benjamin Bustos 


© Springer-Verlag 2011 


Abstract With the increasing amount of 3D data and the 
ability of capture devices to produce low-cost multimedia 
data, the capability to select relevant information has be- 
come an interesting research field. In 3D objects, the aim 
is to detect a few salient structures which can be used, in- 
stead of the whole object, for applications like object regis- 
tration, retrieval, and mesh simplification. In this paper, we 
present an interest points detector for 3D objects based on 
Harris operator, which has been used with good results in 
computer vision applications. We propose an adaptive tech- 
nique to determine the neighborhood of a vertex, over which 
the Harris response on that vertex is calculated. Our method 
is robust to several transformations, which can be seen in the 
high repeatability values obtained using the SHREC feature 
detection and description benchmark. In addition, we show 
that Harris 3D outperforms the results obtained by recent 
effective techniques such as Heat Kernel Signatures. 


Keywords 3D interest points detection - Local features - 
Harris operator 
1 Introduction 


Many applications have benefited from the wide diffusion of 
3D models. Areas such as medicine, engineering, entertain- 
ment, and so on are increasingly relying in processes that 
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involve this kind of information. Coupled with this, an im- 
proved ability of capture devices has been observed, allow- 
ing to generate low-cost three-dimensional objects and make 
extensive use of them. In addition, tasks involving 3D mod- 
els such as mesh analysis and processing are active research 
areas. 

For the same reasons, there is a growing interest in the 
use of detailed models for better representation. This im- 
plies a higher necessity for feature detection tasks. As with 
images, the better is the resolution of a 3D object, the better 
the representation of some entity and therefore it is neces- 
sary to be able to select distinctive points on a 3D model 
in order to keep the efficiency in the processes applied on 
them. Some tasks that benefit from this capability are object 
registration [5], object retrieval and matching [9], mesh sim- 
plification, viewpoint selection [14], and mesh segmentation 
[11, 25], just to name a few. 

The interest point detection on 3D data is a challeng- 
ing problem for several reasons. First, there is no consen- 
sus about the definition of an interest point. A commonly 
used definition (that we use in this paper) relates the measure 
of interest with the level of protrusion of outstanding local 
structures. So, vertices on smooth or nearly planar sections 
of a surface will have low interest, as opposite to vertices 
in regions with uncommon local structure. For instance, in 
a human-shaped model, an interest point detector should se- 
lect vertices on the face, hands, and feet. Second, the topol- 
ogy in 3D meshes is arbitrary. That is, a vertex can have 
an arbitrary number of neighboring vertices. This makes 
the tasks of selecting a local neighborhood around a vertex 
harder. In addition, this drawback causes that different tes- 
sellations can represent the same locality and therefore, an 
interest point detector should be able to deal with that. Third, 
without a well-defined topological structure for meshes, the 
extent of a locality in which a vertex is an interest point is 
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Fig. 1 This figure shows the 
steps we propose to detect 
interest points in 3D meshes 
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unknown or difficult to compute. Finally, there is no addi- 
tional information other than the position of vertices and the 
connectivity information among them. This fact complicates 
the process because the level of interest needs to be mea- 
sured using the available information, which also depends 
of the topology of the mesh. 

In a different scenario, interest point detectors for im- 
ages have reached an acceptable effectiveness. The reason is 
that image structure is well-defined and interest points corre- 
spond with pixels that represent interesting structures in the 
scene captured in the image. In that sense, several methods 
have been proposed for detecting interest points taking into 
account different scales and transformations. As a result, the 
range of applications in computer vision that makes use of 
interest points detection has increased considerably in recent 
times such as image matching [17], image stitching [2], and 
human activity recognition [12, 13], just to name a few. 

However, despite the success of these techniques in the 
image domain, trying to adapt them to 3D meshes is not a 
trivial task. Firstly, the adaptation is not direct because 3D 
meshes structure is very different from images. Secondly, 
the transformations required to be robust in 3D domains 
(isometry, topology, change of tessellations, downsampling, 
etc.) are also different. 

In this paper, we present an effective and efficient exten- 
sion of the Harris operator for 3D meshes. We chose the Har- 
ris operator for several reasons. First, the computation of the 
operator is an efficient and simple task. This is an important 
issue if we want use the interest point detection as a pre- 
liminary stage of subsequent process such as shape match- 
ing, registration or object recognition. Second, Harris-based 
methods have been effectively used in a number of applica- 
tions and they have a high effectiveness as reported in the 
evaluation reports [19, 21]. Finally, an interesting evidence 
has been found recently, which greatly favors the Harris in- 
terest point detection method. Loog and Lauze [16] recently 
showed an important connection between the Harris oper- 
ator and the computational visual attention model of visual 
perception. Roughly speaking, their model proved that inter- 
est points detected by the Harris method have low probabil- 
ity of appearing in other locations in the same image. These 
reasons encouraged us to investigate an effective extension 
of the Harris operator for 3D meshes, trying to comply with 
the robustness to the transformations in 3D domain. 

Our method is outlined in Fig. 1. Given a 3D mesh, not 
necessarily manifold, the main process is performed in a 
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vertex-wise manner. The overall process consists of four 
steps. Firstly, our algorithm determines a local neighbor- 
hood around a vertex. The subsequent tasks are performed 
over this local neighborhood. Secondly, the neighborhood is 
processed so that it is prepared for a fitting step. We try to fit 
a quadratic surface to the set of points. This surface is a good 
representation of the locality and we consider it as an local 
image. Thirdly, we propose to calculate derivatives using a 
smoothing over the surface. We can use these derivatives for 
computing the Harris response for each vertex. Finally, our 
method selects the final set of interest points. A preliminary 
version of our method has been presented in a conference 
paper [23]. 

The contributions of this paper are summarized as fol- 
lows: 


e We improve the process for calculating the Harris opera- 
tor for 3D meshes, making it robust to noise, change of 
tessellations and other transformations which deform the 
mesh structure such as local scaling, shot noise, presence 
of holes, just to name a few. 

e We propose a novel method to define the neighborhood 
size of a vertex, depending on its surrounding structure. 

e We give several options to select a few interest points us- 
ing the information that the Harris operator provides. 

e We present a comprehensive experimentation, trying to 
investigate the effect of different parameters choices and 
comparing our method with effective methods in the state 
of the art. 


The organization of this paper is as follows. Section 2 
presents the related works. Section 3 presents a detailed de- 
scription of our method. Section 4 presents and discusses the 
experimental results. Finally, Sect. 5 concludes the paper. 


2 Related work 


The interest point detection topic emerged in the computer 
vision community with the aim of reducing the amount of 
information used in high-level vision tasks. A pioneering 
work was presented by Harris and Stephens [7], which was 
the basis for many later works. For readers interested in in- 
terest points detectors on images, we recommend the evalu- 
ation paper presented by Schmid et al. [21], which contains 
detailed descriptions and performance evaluation of several 
proposed methods. 
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For 3D meshes, several approaches have been proposed, 
most of which have tried to extend the detectors proposed 
for images. After the SIFT method proposed by Lowe [17], 
a number of extensions have been presented which use 
Difference-of-Gaussians (DoG) as interest point detector. 
Castellani et al. [3] applied the DoG detector over vertices 
in scale space obtained with successive decimations of the 
original shape. Vertices with high response in their DoG op- 
erator are selected as interest points. In the same way, Zou 
et al. [27] proposed to build a geodesic scale space, and sub- 
sequently to apply DoG detector on that space for detect- 
ing interest points on a surface. Also, Zaharescu et al. [26] 
assumed that the vertices of an 3D object have associated 
information such as curvature or photometric properties. 
Defining a discrete Difference-of-Gaussians operator, the 
authors applied this operator on the function defined by 
the associated information over a manifold. This approach 
showed good results in matching of 3D models sequences. 

On the other hand, the geometric diffusion theory can be 
used for detecting interest points on surfaces. The diffusion 
process reveals interesting characteristics from the intrinsic 
geometry of a surface which can be exploited to detect out- 
standing structures. As a 3D surface property related to the 
diffusion process on a manifold, the Laplace—Beltrami op- 
erator has been also used to detect interest points. Hu and 
Hua [9] defined the geometric energy of a vertex as function 
of the eigenvalues and eigenvectors of the Laplace—Beltrami 
spectrum of a given object. Vertices where the energy is a 
maximum are considered as interest points. In addition, the 
energy provides the scale where the selected vertices are in- 
teresting. The selected interest points were used in a match- 
ing task with promising results. On the other hand, Sun et al. 
[24] defined the Heat Kernel Signature as a temporal domain 
restriction of the Heat Kernel on a manifold, which is re- 
lated to the Laplace—Beltrami spectrum. In 3D meshes, each 
vertex has an associated signature. A vertex is selected as 
interest point, when for large time values, its signature has a 
maximum with respect to the neighbor vertices. 

Similarly, Zou et al. [28] proposed to build a scale space 
of the surface geometry via the surface Ricci flow which 
satisfies a set of desired properties for a multi-scale repre- 
sentation. The authors applied the Ricci flow over a metric 
of the surface based on edge lengths. The scale space is rep- 
resented as a matrix with curvature values calculated from 
the set of diffused metrics. Then, the Laplacian of a vertex 
is computed using the cotangent schema using the curvature 
as associated values. A vertex is considered as an interest 
point if it is an extreme of the Laplacian in the 1-ring neigh- 
borhood and neighboring scales. 

Also, curvature-based methods have been proposed. 
Gelfand et al. [5] described an interest point detector based 
on a new descriptor called the integral volume descriptor. 
For each vertex, the amount of volume in the intersection 


of a ball centered in the vertex and the 3D object describes 
an interesting local measure. The authors showed that this 
quantity is closely related to the curvature in the vertex. Ver- 
tices with uncommon integral values are selected as interest 
points. Also using curvature in vertices, Ho and Gibbins [8] 
suggested a measure called the curvedness measure in or- 
der to describe the geometric information in a vertex. The 
curvedness is calculated from the principal curvatures of a 
vertex. This measure can be calculated in different scales by 
selecting different neighborhood sizes which are used to fit 
quadratic patches over which curvatures are computed. Ver- 
tices with extremal values in its curvedness, with respect to 
neighboring points and scales, are selected as interest points. 

Differently, Liu et al. [15] proposed a Monte Carlo strat- 
egy to select a random set of points on a surface with each 
point having the same probability to be chosen. These points 
were used in partial shape retrieval. The assumption be- 
hind this proposal is that the vertices of a shape are sam- 
ples of the original surface and the tasks that use them can 
be affected by shape tessellations. Similarly, Shilane and 
Funkhouser [22] considered random points on a 3D surface, 
selecting only those points that contribute to improve the re- 
trieval performance. With a training phase, it was possible to 
assign a predicted distinction value to each selected point in 
the 3D collection and thus, using that values to assign new 
ones to points of a new shape. 

As another approach, the mesh saliency defined by Lee 
et al. [14] has proven to be a robust feature to many 3D ap- 
plications. The process to compute the mesh saliency of a 
3D object begins calculating a Gaussian-weighted average 
of the mean curvature on a surface. Each vertex in an ob- 
ject is thus associated with the difference of such average 
in different scales, which is the saliency of that vertex. Ver- 
tex with the highest saliency can be considered as interest 
points. 

Conformal parameterization has also been used to pro- 
pose interest points detectors. Methods based on conformal 
parameterization [10, 20] transform a 3D surface into a 2D 
parameterization that can be seen as an image. A 3D to 2D 
mapping is said to be conformal if angles are preserved. 
Once an image is computed, interest points can be detected 
on it, and subsequently these are mapped back to the 3D 
domain. 

Finally, Mian et al. [18] related the repeatability of key- 
points (extracted from partial views of an object) with a 
quality measure based upon principal curvatures. 


3 Interest points detection 
Harris and Stephens [7] proposed an interest points detector 


for images. Their method is a popular technique due to its 
strong invariance to rotation, scale, illumination variation, 
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and image noise [21]. The Harris detector is based on the lo- 
cal autocorrelation function of a signal, which measures the 
local changes of the signal with patches shifted by a small 
amount in different directions. The local autocorrelation is 
defined as: 


ex, y) = Ý WO, vi) [Fi + Ax, yi + Ay) Ey 


Xi, yi 
a) 


where I (-,-) denotes the image function and (x;, y;) are the 
points in the Gaussian function W centered on (x, y), which 
defines the neighborhood area in analysis. 

Using a Taylor expansion truncated to the first order 
terms to approximate the shifted image, we obtain: 


2 
ek Dry Wi; Diy W.Iy D ST 
en W.I,.1y ae Wit; (2) 
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where S = [Ax Ay] is a shift vector, J, and 7y denote the 
partial derivatives in x and y, and along with W are evalu- 
ated in (x;, yj) points. 

Harris and Stephens proposed to analyze the eigenvalues 
of matrix E, which contains enough local information re- 
lated to the neighborhood structure. In addition, to avoid the 
expensive eigenvalue calculation, they proposed to assign to 
each pixel in the image the following value: 


h(x, y) = det(E) — k(tr(E))” (3) 


with k constant. 

The Harris operator has been used in many applications 
in image processing and computer vision by its simplicity 
and efficiency. However, the problem with 3D data is that 
the topology is arbitrary and it is not clear how to calculate 
the derivatives. To cope this problem, Glomb [6] suggested 
some approaches. We take this work as a basis for proposing 
a robust interest points detector on 3D meshes. 


3.1 Robust Harris operator on 3D meshes 


Given a vertex of a 3D object, we are interested in calcu- 
lating the Harris operator value associated with that point. 
A 3D object is represented as a set of vertices V and a set of 
faces F with adjacency information between these entities. 
In addition, our method is not restricted to manifold meshes. 

Let v be the analyzed vertex and V;.(v) the neighborhood 
considering k rings around v. Figure 2 shows vertex v (black 
circle), the first ring around v (path formed by green circles), 
the second ring (path formed by blue circles), and the kth 
ring (path formed by yellow circles). All these vertices cor- 
respond to the neighborhood V;(v). The method to calculate 
k will be explained later in this section. 
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Fig. 2 Point v and its neighbor rings. Firstly, Vı (v) is composed by 
vertices connected by strong edges. Secondly, V2(v) is composed by 
vertices up to those connected by dashed edges. Finally, Vg (v) is com- 
posed by all vertices until those connected by pointed edges 


We calculate the centroid of Vg(v) and translate the set 
of points so the centroid is in the origin of the 3D coordi- 
nate system. Then, we compute the best fitting plane to the 
translated points. To do so, we apply Principal Component 
Analysis to the set of points and we choose the eigenvec- 
tor with the lowest associated eigenvalue as the normal of 
the fitting plane. In our opinion, applying PCA is a better 
choice than the least square fitting because the assumption 
z= f (x, y) does not have a good behavior when the data do 
not exhibit such functional characteristic. 

The set of points is rotated so that the normal of the fit- 
ting plane is the z-axis. As we choose the less principal 
component as normal, the points exhibit a good spread in 
theX Y-plane after rotation and therefore we can only work 
in XY-plane to calculate the derivatives. As final step before 
calculating derivatives, we translate the set of points so that 
the point v is in the origin of the XY-plane. This step will 
facilitate the further analysis. 

To calculate derivatives, we fit a quadratic surface to the 
set of transformed points. Using least square method, we 
find a paraboloid of the form: 


z= f(x, y)= Bx + pary + Sy? + pax + psy + po: (4) 

We chose a quadratic surface with only six terms because 
it represents a paraboloid. That is, it is the best choice if 
we need a function of two variables with quadratic terms. 
Adding more terms implies that it is possible to fit a more 
complex surface. However, more complex surfaces do not 
have well-defined derivatives in certain points of the do- 
main. In addition, we needed a simple expression in order 
to apply the derivatives. 

AS we are interested in derivatives in the point v, one 
could directly evaluate the derivatives of f(x, y) in the point 
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The above expressions should be a good estimate of 
derivatives. However, these can be influenced by noise. In- 
stead, we propose to apply a Gaussian function as proposed 
originally by Harris and Stephens [7]. However, a difficulty 
arises because in the original expression the derivatives are 
discrete functions and our derivatives are continuous func- 
tions. To address this problem, we propose to apply the inte- 
gration of the derivatives with a continuous Gaussian func- 
tion as follows: 
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where ø is a constant, which defines the support of the 
Gaussian function and the factor 1/v 27o is a normaliza- 
tion value. 
Using calculus, we can reduce the expressions to 


C 


A= py +2p +2p2, a0) 
B = ps +2p3 +2p3, (11) 
C = paps + 2pı p2 + 2p2 p3. (12) 


Finally, we can formulate the matrix E associated with 
the point v using the previously calculated values: 


A C 
r (13) 


The Harris operator value in the point v is calculated as 
in (3). 


3.2 Adaptive neighborhood size 


Several approaches can be considered to select the number 
of rings around a point as neighborhood. If the object tes- 
sellation is uniform, i.e., almost all triangles in the object 
have the same size, we can use a constant number of rings 
to all points, or use the points contained in a ball of radius r 
and centered in point v. However, in irregular and complex 
meshes, these methods do not approximate a neighborhood 
adequately. 


To tackle this problem, we propose an adaptive tech- 
nique. Our method selects a different neighborhood size de- 
pending on the tessellation around a point. Let us consider 
an object as a graph G(V’, E’), where V’ = V and E’ is the 
set of edges obtained from the adjacency information of the 
object. 

Given a point v € V’, a k-ring around v is the set of points 
where the length of the shortest path to v is k: 


ring; (v) = fw eV’ | shortest_path(v, w)| = k}. (14) 
The distance from a point v to the ring; (v) is defined as 


max 
wering; (v) 


dring(v, ring, (v)) = |v — wllo. (15) 


Finally, we define the neighborhood size of a point v as 
radius, = {k EN, dring (v, ring, (v)) > 6 and 
dring(v, ring,_(v)) < ô}, (16) 


where ô is a fraction of the diagonal of the object bounding 
rectangle. 

It is important to note that the proposed method always 
finds a neighborhood to a point, even with complex and ir- 
regular tessellations around that point. 


3.3 Selecting interest points 


With each vertex associated with its Harris operator value, 
we propose two ways to select the interest points of a given 
object. Firstly, we preserve the vertices which are local max- 
imum. To do so, we select a vertex v which holds the follow- 
ing condition: 


h(v)>h(w), Vw € ring; (v). (17) 


Secondly, we propose two approaches to select the final 
set of interest points. 


e Select the points with the highest Harris response. We can 
pick a constant fraction of interest points depending on 
the application. In this proposal, we obtain the points with 
higher saliency and therefore, some portions of the object 
do not have interest points. 

e Representatives of Interest Points Clusters. This approach 
can be used when we want a good distribution of inter- 
est points in the object surface. This proposal consists of 
two steps. First, we sort the pre-selected interest points 
according to their Harris operator value in decreasing or- 
der. Second, we apply Algorithm 1 to cluster the sorted 
points and select the final set of interest points. 

The value of p can be considered as a fraction of the 
diagonal of the object bounding rectangle and it has effect 
in the number of returned interest points. 


Figure 3 shows the result of the two options to select in- 
terest points. 
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Fig. 3 Selection options. (a) Armadillo model. (b) Selected points with highest Harris response. (c) Selected points by clustering 


Algorithm 1 Interest Points Clustering 

Require: Set P of pre-selected interest points in decreasing 
order of Harris operator value 

Ensure: Final set of interest points 


1: Let Q be a set of points 

2 O<G 

3: fori < 1 to |P| do 

4: if minjeri,joy || Pi — Qj ll2 > p then 
5: Q < QU{P;} 

6: end if 

7: end for 

8: Return Q 


4 Experimental evaluation and discussion 


In this section, we show the experimental results of our 
implementation of the Harris 3D method using a standard 
benchmark. The presentation of results is divided in two 
parts. First, we experiment with several aspects and parame- 
ters of our method. The objective is to investigate the effect 
of the parameters on the repeatability of interest points. Sec- 
ond, we compare our method with methods in the state of 
the art, namely the Heat Kernel Signatures proposed by Sun 
et al. [24] and the Salient Points detection method proposed 
by Castellani et al. [3]. 


4.1 The data set 


For our experiments, we use the SHREC robust feature de- 
tection and description benchmark [1], which is available in 
the public domain. This data set is composed of triangular 
meshes with approximately 10,000—50,000 vertices. 

The collection consists of three basic shapes (null shapes) 
from which a set of transformations have been applied. For 
each null shape, nine transformations were used: isome- 
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try (non-rigid transformation), topology, big holes and mi- 
cro holes, local and global scaling, additive Gaussian noise, 
shot noise, and downsampling (less than 20% of the original 
points). 

Each transformation is performed in five versions. In all 
transformations, except isometry and scale, each version 
represents a strength level, so high levels correspond with 
stronger transformations. The scales used in scaling trans- 
formation were 0.5, 0.875, 1.25, 1.625, and 2. For the isom- 
etry transformation, each version reflects a non-rigid trans- 
formation of the null shape. Therefore, each null shape has 
45 transformed shapes and there are 138 shapes in the whole 
collection. 


4.2 Evaluation methodology 


An interest point detection method returns a set of detected 
points F(Y) = {yx}, for each shape Y (typically, |F(Y)|<« 
|Y |). The performance is measured by comparing the inter- 
est points computed for transformed shapes and the corre- 
sponding null shape. 

The quality of the interest points detection was measured 
using the repeatability criterion. Assuming for each trans- 
formed shape Y in the data set the ground-truth dense cor- 
respondence to the null shape X to be given in the form 
of pairs of points Co(X, Y) = {y}, xx)} "an (and same way, 
Co(Y, X)), an interest point yg € F(Y) is said to be repeat- 
able if a geodesic ball of radius R around the correspond- 
ing point x; : (x, yk) € Co(X, Y) contains an interest point 
xj € F(X). The subset F, (Y) € F(Y) of repeatable interest 
points is given by 


Fr x (Y) = {yk € F(Y) : F(X) N Br) 49, 
(xp, yk) € Co(X, Y)}, (18) 


where Br(x;,) = {x € X : geod(x,x,) < R} and geod de- 
notes the geodesic distance function in X. The repeatability 
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Fig. 4 (a) Repeatability values obtained by selecting spatial neighbor- the range [0.01, 0.1] for the radii. (b) Average repeatability for spatial 


hoods. The radii of the ball was chosen as a fraction of the diagonal of neighborhoods with radii 0.02 and 0.03, respectively. The comparison 
the bounding box of the object. Results show average repeatability for was done by levels of transformation 

rep(Y, X) of F(Y) in X is defined as the percentage of in- It is important to point out that for all experiments we 

terest points from F(Y) that are repeatable: used a basic configuration of parameters and only the ana- 

lyzed parameter was changed. Our basic configuration con- 

rep(Y, X) = |Fr.x (Y)| f (19) sisted of adaptive local neighborhood with ô = 0.01, K = 

|F(Y)| 0.04, and selection of 1% of the number of vertices with 


higher Harris response as final interest points set. All graphs 
presented in our results were adequately scaled for a better 
visualization. 


For a transformed shape Y and the corresponding null 
shape X, the overall feature interest points detection quality 
was measured as (rep(Y, X) + rep(X, Y))/2. The value of 
R = 5 was used in the benchmark. This radius constitutes 
approximately 1% of the shapes diameter. Interest points 
without ground-truth correspondence (e.g. in regions in the 
null shape corresponding to holes in the transformed shape) 
were ignored. 


4.3.1 Local neighborhood 


In order to asses the importance of the determination of 
the local neighborhood around an analyzed vertex, we per- 
formed three experiments: spatial, adaptive and ring neigh- 
4.3 Analysis of parameter values Hook: 
Spatial neighborhoods Our first experiment consisted in 
evaluating the repeatability of our method when local neigh- 
borhoods were determined as the set of vertices lying inside 
of a ball centered in the analyzed vertex. We varied the ra- 
dius of the ball and calculated the average repeatability for 
e The type and size of local neighborhood. We tested three all transformations and for all levels. The radii were taken as 
options: spatial neighborhood, adaptive neighborhood fraction of the diagonal of the bounding box of the object. 
and ring neighborhood. We are interested in evaluating Figure 4a shows the results of this experiment. 


In this section, we present the experimental results of our 
method. We experimented with several aspects and investi- 
gated the effect of the parameters on the repeatability of our 
proposal. Specifically, we evaluated the following aspects: 


the effectiveness of each option. We varied the fraction of the diagonal in the range 
e The parameter K. We tested with different values for this  [0.01, 0.1]. The highest values of repeatability were ob- 
parameter in order to investigate its effect in the calcula- tained for 0.02 and 0.03. In order to find the best value, we 
tion of Harris response. plot the repeatability curves for these values with respect to 
e The type of interest point selection method. We only eval- the level of transformation (see Fig. 4b). Both curves are 
uated the method that selects the points with higher re- very similar, with a slight advantage for fraction = 0.03. 
sponse. We intended to figure out the effect of the number There are two important issues that we should remark. 
of selected interest points in the effectiveness of our algo- First, the average repeatability decreases with large neigh- 
rithm. borhoods. This is because large neighborhoods cannot be 


a Springer 


I. Sipiran, B. Bustos 


Fig. 5 Average repeatability for 100 T 
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Table 1 Repeatability of our method using spatial neighborhood with 
fraction = 0.03. Average number of detected points: 303 


0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 
Size adaptive neighborhood 


Table 2 Repeatability of our method using adaptive neighborhoods 
with ô = 0.01. Average number of detected points: 303 


Transform. Strength 
1 <2 <3 <4 <5 

Isometry 75.77 79.89 78.01 79.45 79.85 
Topology 76.25 76.39 76.21 76.19 76.18 
Holes 75.56 75.52 75.36 75.34 74.98 
Micro holes 76.19 76.17 76.10 75.98 75.92 
Scale 82.58 80.01 77.38 74.79 72.29 
Local scale 74.33 72.49 69.34 65.46 62.66 
Sampling 73.21 71.06 69.35 66.75 59.51 
Noise 73.81 73.14 71.78 69.94 67.48 
Shot noise 76.32 76.36 76.04 76.09 75.90 
Average 76.00 75.67 74.40 73.33 71.64 


Transform. Strength 
1 <2 <3 <4 <5 

Isometry 93.59 95.14 94.43 95.05 95.18 
Topology 93.48 93.43 93.25 93.16 93.15 
Holes 92.08 92.04 91.60 91.25 90.67 
Micro holes 93.59 93.59 93.56 93.55 93.51 
Scale 94.29 93.80 93.42 93.04 92.60 
Local scale 93.49 92.89 91.27 88.81 86.55 
Sampling 92.23 90.40 87.83 84.52 77.98 
Noise 92.33 81.83 72.51 66.75 63.44 
Shot noise 93.54 92.60 91.27 89.93 88.20 
Average 93.18 91.75 89.90 88.45 86.81 


well fitted by a quadratic surface, and therefore the calcu- 
lation of derivatives and the Harris response are not robust. 
Second, as it can be seen in Fig. 4b, repeatability always 
decreases with stronger levels of transformations (except in 
those where the level does not represent the strength of the 
transformation). Nevertheless, this behavior was expected. 
Table 1 shows the repeatability for each transformation 
and each level for fraction = 0.03. Although there are trans- 
formations where, on average, the repeatability is greater 
than 75% (for example, isometry, topology, holes, micro 
holes, scale and shot noise), spatial neighborhoods have 
some problems. On the one hand, vertices belonging to a 
spatial neighborhood do not necessarily belong to the same 
local surface around the analyzed vertex. Moreover, the set 
of vertices could not belong to the same connected compo- 
nent. Having neighborhoods with vertices not belonging to 
the local surface, the fitting task is not robust, damaging the 
overall process and therefore the effectiveness. On the other 
hand, spatial neighborhoods do not ensure a good balance 
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around the analyzed vertex, which is a problem if the mod- 
els have poor triangulations. 


Adaptive neighborhoods Our second experiment consisted 
in varying the parameter ô with the adaptive local neigh- 
borhood approach. Figure 5 shows the average repeatability 
obtained in this experiment. Table 2 shows all values of re- 
peatability obtained by the best configuration (ô = 0.01). We 
tested values smaller than 0.01 for ô; however, the improve- 
ment was not significant. 

Similarly to the spatial neighborhoods, the average re- 
peatability decreases as the extent of the neighborhoods in- 
creases. Large neighborhoods cause an inadequate fitting, 
affecting the effectiveness of the method. Results show an 
improvement of more than 15% in almost all entries with 
respect to those obtained by spatial neighborhoods. This 
percentage represents approximately 45 repeatable interest 
points detected in difference, compared to the spatial neigh- 
borhood effectiveness. However, although the results im- 
prove significantly, there is still visible a rapid fall in the re- 
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peatability values for noise transformation, even below those 
obtained by spatial neighborhoods. 

From Table 2 we can note some interesting issues. On 
the one hand, from the transformations where the level rep- 
resents strength, some of them present only a slight fall of 
average repeatability between the minimum and maximum 
levels: topology (0.33%), holes (1.41%) and micro holes 
(0.08%). On the other hand, from stronger level, the worst 
performance is obtained by sampling and noise transforma- 
tions. 

The behavior with the sampling transformation is ex- 
pected because our method relies on selecting a fraction 
of the number of vertices as interest points. So the num- 
ber of interest points detected in level 1 is much larger than 
in level 5. Clearly, the number of repeatable interest points 
in this transformation cannot be larger than the number of 
interest points in high levels. Therefore, the repeatability 
decreases considerably with downsampling. Also, with re- 
spect to the noise, the level of distortion of the meshes in 
stronger levels of this transformation is high, causing a con- 
siderable deformation in the shapes. Surely, local neighbor- 
hoods are affected considerably and the process of fitting is 
not robust. We argue that the decrement with respect to spa- 
tial neighborhoods is due to the use of geodesic distances 
when collecting the rings around the analyzed vertex. As 
the noise is applied in the direction of the normal of each 
vertex, the same parameter ô in different levels of transfor- 
mations can determine neighborhoods of considerable dif- 
ferent sizes. Obviously, for the reasons explained before, the 
repeatability of our method is affected. 


Ring neighborhoods The third experiment evaluated the 
repeatability with different numbers of rings selected as 
neighborhood. Figure 6 shows the result for this experiment. 
Table 3 shows all repeatability values. 

From Fig. 6 we can observe the same effect on the size 
of the neighborhood. With larger neighborhoods, average 
repeatability begins to decrease systematically. The inter- 
esting thing about this experiment is to note that a very 
small neighborhood gave the best results. One reason could 
be that, with small neighborhoods, the fitting step is ro- 
bust and, thus, derivatives are well calculated. Nevertheless, 
small neighborhoods are unstable in presence of noise or an- 
other distortion transformation. Still, repeatability values for 
stronger levels of noise improve up to 13% with respect to 
spatial neighborhoods and up to 15% with respect to adap- 
tive neighborhoods. 


4.3.2 Parameter K 
The parameter K is used in (3) to calculate the Harris 


response for a given vertex. This parameter needs to be 
tuned experimentally. We varied the parameter in the range 


Table 3 Repeatability of our algorithm using one ring neighborhood. 
Average number of detected points: 303 


Transform. Strength 
1 <2 <3 <4 <5 

Isometry 95.86 96.74 96.51 96.78 96.79 
Topology 95.92 95.92 95.84 95.75 95.70 
Holes 93.60 93.77 93.74 93.83 93.48 
Micro holes 95.86 95.86 95.88 95.89 95.93 
Scale 96.76 96.54 96.24 95.93 95.39 
Local scale 95.92 94.97 93.49 91.35 89.01 
Sampling 94.69 93.55 92.19 89.05 80.26 
Noise 92.97 92.07 90.82 89.71 88.54 
Shot noise 95.99 95.57 94.90 93.62 92.41 
Average 95.28 95.00 94.40 93.55 91.95 


[0.04, 0.1] and tested the average repeatability. Figure 7 
shows the results. In our implementation, the best result 
was for K = 0.07. However, the improvement was not sig- 
nificant with respect to our default value K = 0.04 (ap- 
proximately 0.34%, representing almost 2 repeatable inter- 
est points of difference). 


4.3.3 Interest point selection 


In Sect. 3.3, we proposed two options for the final selec- 
tion of interest points. The clustering approach is interesting 
from the point of view of applications requiring points well 
distributed over the whole surface. Nevertheless, as the pro- 
cess is based on grouping vertices with high responses, this 
step is not necessarily robust according to the repeatability 
criterion. Therefore, we did not consider this approach in 
our evaluation. 

On the other hand, selection of vertices with higher re- 
sponse is interesting because the number of selected vertices 
is an important issue in applications. For example, in shape 
matching, we might be interested in only a few points, as 
the efficiency of matching is closely related to the number 
of points to be matched. So, we evaluated the effect of re- 
ducing the number of selected vertices. 

Our method selects a fraction of the number of vertices 
as interest points, so the smaller is the fraction, the fewer in- 
terest points are selected. We varied the fraction in the range 
[0.001, 0.01]. Figure 8 shows our results. 

Clearly, the average repeatability increases as the number 
of interest points is increased. This trend is maintained for 
values of fractions larger than 0.01. An important aspect to 
note is that by reducing the number of interest points in half, 
with respect to the value of fraction 0.01, the average re- 
peatability decreases approximately in 6.34%. This is an ad- 
vantage because if we need to apply subsequent processes, 
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Fig. 6 Average repeatability 
with ring neighborhoods sizes 
taken from the range [1, 10] 
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we can opt for reducing the number of interest points se- As an additional test, we combine the parameter values 


lected, at expense of slightly reducing the robustness of the with the best effectiveness, obtaining a slight improvement. 
delivered points. The number of interest points will finally | Table 4 shows the repeatability values for all transforma- 
depend on the application and the trade-off between robust- tions in all levels. Compared to Table 3, on average results 
ness and efficiency. improve for all levels. 
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Table 4 Repeatability of our algorithm using combination of values 
with the best effectiveness. Average number of detected points: 303 


Table 5 Repeatability of HKS1 feature detection algorithm. Average 
number of detected points: 35 


Transform. Strength 
1 <2 <3 <4 <5 

Isometry 96.01 96.73 96.26 96.60 96.62 
Topology 96.01 95.97 95.82 95.73 95.71 
Holes 94.62 94.43 94.10 94.01 93.81 
Micro holes 96.01 96.01 95.98 95.96 95.95 
Scale 97.06 96.89 96.28 95.62 94.94 
Local scale 96.24 94.96 93.40 91.26 88.84 
Sampling 95.31 93.62 92.08 89.13 80.42 
Noise 93.09 92.58 91.59 90.33 88.79 
Shot noise 96.03 95.66 95.00 93.83 92.79 
Average 95.60 95.21 94.50 93.61 91.99 


Transform. Strength 
1 <2 <3 <4 <5 

Isometry 83.20 87.39 88.44 87.29 87.21 
Topology 79.73 81.27 81.49 81.08 80.94 
Holes 80.15 77.57 75.86 73.35 71.16 
Micro holes 81.02 80.76 80.68 80.43 80.50 
Scale 79.99 79.78 79.90 80.18 80.36 
Local scale 80.38 80.65 78.84 75.55 72.99 
Sampling 82.70 82.23 82.66 82.04 78.99 
Noise 75.80 74.55 72.37 69.34 68.23 
Shot noise 79.96 81.14 81.15 80.07 78.77 
Average 80.33 80.59 80.15 78.82 771.68 


4.4 Comparison with other methods 


In order to compare our method with the state of the art, 
we selected two recent methods for detecting interest points 
on 3D meshes: Heat Kernel Signatures [24] and Salient 
Points [3]. Next, we specify the configuration used for these 
methods: 


Heat kernel signature (HKS) As this method relies on the 
eigendecomposition of the Laplace—Beltrami operator on 
the mesh, we needed to simplify the meshes to approxi- 
mately 10,000 vertices (we used Garland’s method [4]). In- 
terest points were computed on the simplified meshes and 
mapped back to the original mesh by using the nearest 
neighbor vertex. We used the value t = 0.1 from the total 
area of the surface to evaluate the Heat Kernel Signature 
and 2-ring neighborhood in order to detect interest points. 
We present three variations: 


e HKS1: No pre-processing step was used. We imple- 
mented this variant based on the original HKS implemen- 
tation.! 

e HKS2: The Geomagic sotfware was used for removing 
non-manifold edges, and faces were consistently oriented. 
Results for this variant were taken from the original re- 
port [1]. 

e HKS3: Filtering using persistent homology was used to 
discard unstable feature points. Results were also taken 
from [1]. 


Salient points (SP) We chose the best performance of this 

method from the SHREC feature detection benchmark. 
Tables 5, 6 and 7 show the repeatability obtained by the 

variants of the Heat Kernel Signature method. Table 8 shows 


‘http://www.geomtop.org/sunjian/software/hks.html. 


Table 6 Repeatability of HKS2 feature detection algorithm. Average 
number of detected points: 23 


Transform. Strength 
1 <2 <3 <4 <5 

Isometry 98.08 98.72 98.01 97.88 98.04 
Topology 97.44 96.10 92.26 91.22 88.64 
Holes 91.48 90.60 86.78 83.73 81.86 
Micro holes 98.08 96.69 96.00 95.52 94.87 
Scale 99.36 99.36 98.50 97.90 97.68 
Local scale 98.08 94.83 90.09 83.05 78.31 
Sampling 97.05 97.88 97.39 96.27 92.35 
Noise 95.30 92.78 91.67 89.24 87.62 
Shot noise 98.08 96.22 93.39 90.45 87.32 
Average 96.99 95.91 93.79 91.70 89.63 


the results of the Salient Points method. All comparison 
were done with respect to our best results shown in Table 4. 
Our method is represented by H3D. 

As can be seen, our method outperformed the HKS1 vari- 
ant without pre-processing and the Salient Points method. 
The benefit is consistent in each entry of the table, which is 
an important result regarding the relevance of the techniques 
compared. With respect to HKS2 and HKS3, these variants 
present significant improvements. However, they need well- 
formed shapes in order to work properly, which affects their 
efficiency considerably. 

Table 9 presents the best method for each entry in the re- 
peatability table. It can be seen that Harris 3D method out- 
performs the rest in stronger levels (4—5). The Heat Kernel 
based method is predominant in isometry, scale and sam- 
pling transformations. In the isometry transformation, the 
effectiveness of this method (HKS3 variant) is perfect with 
an average repeatability of 100%. A similar scenario takes 
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Table 7 Repeatability of HKS3 feature detection algorithm. Average 
number of detected points: 23 
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Table 9 Methods with the best performance by transformations and 
strengths. HKS—Heat Kernel Signature. H3D—Harris 3D 


Transform. Strength 
1 <2 <3 <4 <5 

Isometry 100.00 100.00 100.00 100.00 100.00 
Topology 94.44 90.38 87.45 88.70 85.76 
Holes 80.54 79.00 75.25 72.10 69.99 
Micro holes 100.00 100.00 98.15 96.58 95.64 
Scale 100.00 100.00 100.00 98.61 97.78 
Local scale 97.44 96.79 93.02 87.25 82.90 
Sampling 100.00 100.00 100.00 100.00 96.20 
Noise 100.00 95.19 93.16 89.37 85.77 
Shot noise 100.00 95.30 90.03 82.10 74.38 
Average 96.94 95.19 93.01 90.52 87.60 


Transform. Strength 
1 <2 <3 <4 <5 

Isometry HKS3 HKS3 HKS3 HKS3 HKS3 
Topology HKS2 HKS2 H3D H3D H3D 
Holes H3D H3D H3D H3D H3D 
Micro holes HKS3 HKS3 HKS3 HKS3 H3D 
Scale HKS3 HKS3 HKS3 HKS3 HKS3 
Local scale HKS2 HKS3 H3D H3D H3D 
Sampling HKS3 HKS3 HKS3 HKS3 HKS3 
Noise HKS3 HKS3 HKS3 H3D H3D 
Shot noise HKS3 HKS2 H3D H3D H3D 
Average HKS3 HKS3 HKS3 H3D H3D 


Table 8 Repeatability of SP algorithm. Average number of detected 
points: 409 


Transform. Strength 
1 <2 <3 <4 <5 

Isometry 86.17 87.42 87.24 87.76 88.15 
Topology 86.18 85.63 85.58 85.56 85.56 
Holes 85.72 85.10 84.34 83.56 82.58 
Micro holes 68.52 62.27 57.96 54.75 51.99 
Scale 89.80 88.28 86.82 85.14 83.70 
Local scale 85.73 84.97 84.48 83.33 82.12 
Sampling 85.02 83.15 82.21 79.94 77.61 
Noise 87.31 85.43 83.28 81.36 79.40 
Shot noise 85.95 84.42 82.77 81.76 81.23 
Average 84.49 82.96 81.63 80.35 79.15 


place in the scale transformation, where for small scales 
(levels 1-2-3) the average repeatability is 100%. For the 
sampling transformation, the average repeatability is near to 
perfect too. 

The reason for the good performance of the Heat Kernel 
based method in the aforementioned transformations is its 
intrinsic definition. This property allows it to appropriately 
define a characteristic shape based on the spectrum of the 
Laplace—Beltrami operator, and therefore it is robust against 
isometry and rigid transformations such as scaling. On the 
other hand, it is also robust to different samplings of the in- 
put mesh, as the interest points are selected for being points 
with large values of Heat Kernel Signatures in large times. 
That is to say, the interest points selection process is robust 
to different tessellations. 

Differently, in transformations that deform the local 
structure of shapes, Harris 3D obtained the best results. 
Firstly, there is a total predominance of our method with 
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respect to the holes transformation. Secondly, our method 
outperforms the rest in most levels (3—5) for topology, local 
scale and shot noise transformations. Finally, it is interesting 
that our method is also the best in stronger levels for micro 
holes and noise transformations. Averages Table 9 represent 
the majority in each level, so we can observe the predomi- 
nance of Harris 3D in levels 4 and 5. 

Figure 9 shows some examples of interest points over a 
class of shape of the SHREC feature detection benchmark 
using Harris 3D. 


5 Conclusions 


We presented a robust interest point detector for 3D meshes. 
This task is of importance due to the ability of reducing the 
amount of information needed in subsequent processes. Our 
algorithm effectively adapts the well-known Harris corner 
detection for images in order to be used for 3D meshes. Our 
proposal has showed to be effective, obtaining high repeata- 
bility values using the SHREC feature detection and descrip- 
tion benchmark. 

Our method is robust for several reasons. First, the use 
of a Gaussian function to smooth the derivatives surface 
contributes to effectively mitigate the effect of local defor- 
mations introduced by noise, holes, change of tessellations, 
and so on. Second, our proposal of adaptive neighborhoods 
improves considerably the alternative of spatial neighbor- 
hoods. In addition, from the results obtained by using ring 
neighborhoods, our experiments confirm the fact that bal- 
anced neighborhoods favor the overall process. This is be- 
cause, with a balanced set of points, the approximation of 
derivatives in the analyzed point is better. Finally, along with 
the task of final selection of interest points, Harris 3D pro- 
poses several alternatives for its effective use, making an in- 
teresting alternative in applications such as shape matching 
and registration, just to name a few. 
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Fig.9 Shapes with interest 
points. Interest points are 
represented with small red balls 
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Furthermore, our performed experiments suggest that our 
method could be used in extreme conditions with good re- 
sults. This is an important advantage with respect to the state 
of the art, since it allows us to deal with several kinds of 
shapes and expand the spectrum of possible applications in 
the future. 
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